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The phase space for graphene’s minimum conductivity cr m i n is mapped out using Landauer theory 
modified for scattering using Fermi’s Golden Rule, as well as the Non-Equilibrium Green’s Function 
(NEGF) simulation with a Monte Carlo sampling over impurity distributions. The resulting ‘fan 
diagram’ spans the range from ballistic to diffusive over varying aspect ratios ( W/L ), and bears 
several surprises. The device aspect ratio determines how much tunneling (between contacts) is 
allowed and becomes the dominant factor for the evolution of cr min from ballistic to diffusive regime. 

We find an increasing (for W/L > 1) or decreasing ( W/L < 1) trend in <7min vs. impurity density, 
all converging around 12Sq 2 /ir 3 h ~ 4 q 2 /h at the dirty limit. In the diffusive limit, the conductivity 
quasi-saturates due to the precise cancellation between the increase in conducting modes from charge 
puddles vs the reduction in average transmission from scattering at the Dirac Point. In the clean 
ballistic limit, the calculated conductivity of the lowest mode shows a surprising absence of Fabry- 
Perot oscillations, unlike other materials including bilayer graphene. We argue that the lack of 
oscillations even at low temperature is a signature of Klein tunneling. 


Since its discovery in the last decade, single layer 
graphene has catalyzed widespread research 1 stemming 
from its extraordinary material properties. Multiple elec¬ 
tronic, spintronic and opto-electronic applications are 
predicted to arise from the entire class of 2D materi¬ 
als emergent in graphene’s footsteps 2 . Despite intense 
scrutiny, there exist many unresolved issues that con¬ 
tinue to make the material fascinating. Among them is 
the physics of the minimum conductivity, cr m i n around 
the Dirac point, where the density of states is expected 
to vanish. Instead of vanishing accordingly, cr m i n for 
a ballistic sheet with large width to length aspect ra¬ 
tio ( W/L 1) is shown to be a universal constant 
gq = 4g 2 /7rh 3,4 . This arises from the preponderance 
of tunneling through a continuum of subbands with near 
zero bandgaps. In these structures (W » L samples), 
a series of exponentially decaying tunnel transmissions 
adds up to an overall Ohmic term that factors out of 
the ballistic conductance G = aW/L. Measured cr m i n s, 
however, are typically in the range 4 — 12 g 2 /h 5,7 ^ 9 , ex¬ 
cept Ref. 3 , larger than gq. This is surprising given that 
these experiments are mostly on dirty samples where we 
expect the conductivity to be not only non-universal, 
but certainly smaller than the ballistic limit. The in¬ 
crease in cr m i n from gq arises from charged impurities on 
the substrates that create electron and hole puddles and 
contribute states to the charge neutrality point 10 . How¬ 
ever an opposite, decreasing trend of cr m i n vs. impurity 
concentration (ni mp ) was demonstrated theoretically by 
Adam et.al. in Ref. 6 within Boltzmann transport theory, 
as well as experimentally in Ref. 5 . Clearly there are sev¬ 


eral disjointed pieces that have yet to come together to 
provide a complete phase picture of the evolution of cr m i n 
with sample quality. 

In this paper, we use quasi-analytical Landauer equa¬ 
tion as well as numerical NEGF (within the Fisher-Lee 
formulation) 14 to map out the entire phase space of cr m i n 
for varying 7ii mp and W/L (Fig. 1). Our results clearly 
show that the missing link is the total tunneling current 
(a function of W/L), a piece of physics typically ignored 
in semi-classical models. The observed quasi-saturation 
arises due to a trade-off between the number of modes 
and the scattering time r from charge puddles, as we 
move from the ballistic to diffusive regime. The total 
conductivity can be written as 
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where Go = 4 q 2 /h is conductance quantum including 
spin and valley degeneracy, M p and M e are the number 
of propagating and evanescent modes and T is the corre¬ 
sponding mode averaged transmission probability. While 
this equation defines an absolute lower bound on conduc¬ 
tivity at gq = Aq 2 /nh (dashed line in Fig. 1 top), we will 
shortly show that for dirty samples with impurity den¬ 
sity ~ 3 - 5 x 10 12 /cm 2 , it predicts a quasi-saturating 
g min ~ 4 q 2 /h, consistent with experiments (Fig. 1). 
Part of the fan diagram for W L , the decreasing 
trend in cr m i n in Fig. 1 obtained earlier using Boltz¬ 
mann transport equation, arises naturally in our model 
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FIG. 1. (Color online)(a) Fan diagram of quasi-analytical 
<7 min for W — 500 nm with varying W/L (inset shows conduc¬ 
tance G). The ballistic a m i n is exactly at <jq — 4q 2 /nh. The 
two new features are (1) quasi-saturation at high impurity 
density to ~ 128q 2 /7Y 3 h and (2) a flip in curvature between 
aspect ratios, (b) NEGF calculated cr m i n averaged over pud¬ 
dle geometries (inset). The data saturate at ~ 4 q 2 /h in dirty 
graphene. Dotted green curve is experimental data from Ref. 5 
and open circles are theoretical predictions from Ref. 6 


is the background doping). For the opposite ballistic 
limit, wide samples have a conductivity that dips down 
to the quantized value gq to generate the rest of the 
fan diagram. At the same time, narrow ballistic samples 
with limited tunneling show a conductance quantization 
Go that bears a spectacular robustness with temperature 
and a remarkable absence of Fabry Perot (FP) resonance 
even at low temperature. We interpret the absence of FP 
(Fig. 3) as a clear signature of Klein tunneling, where the 
linear relativistic electron transmits perfectly at normal 
incidence due to pseudospin conservation, contrary to the 
prediction of nonrelativistic Schrodinger equation (which 
applies to Bilayer graphene as we show). Our results are 
supported by numerical NEGF sampled over a Monte 
Carlo distribution of charged impurities. 

Modeling charged impurities: The lack of dangling 
bonds makes direct chemisorption of charged impuri¬ 



FIG. 2. Averaging (a) the pristine graphene density of states 
with (b) a normal distribution of random potentials (c) erases 
the Dirac point, (d) The variance of the Gaussian is calcu¬ 
lated self-consistently and refitted with a simplified expres¬ 
sion, closely matching with the self-consistent calculation in 
the dirty limit. 


ties difficult on graphene. However, dielectric substrates 
can have charged impurities that play a significant role 
on transport around graphene’s Dirac point. The ph- 
ysisorption of charged impurities randomly dopes the 
graphene, creating a Gaussian distribution in energy of 
Dirac points around neutrality. The resulting erasure of 
the Dirac point is already seen in quantum capacitance 
measurements 11 . We can average the linear density of 
states of graphene (Fig. 2a) over a Gaussian distribution 
of potentials (Fig. 2b), with average potential zero, vari¬ 
ance ct£, and potential Ui at the ith location. The exact 
expression involving error functions was worked out by 
Li et. al . 12 , but we can express it in a simpler form that 
interpolates between the low-energy parabolic and high 
energy linear behavior (Fig. 2c). 
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Eq. 2 shows that the variance ge has a direct impact on 
the minimum density of states. Fig. 2d shows that ge in¬ 
creases with charged impurity concentration, so that the 
minimum number of modes for conduction is proportional 
to the statistical variance of charge impurities. This has 
also been worked out by solving Poisson’s equation in 
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cylindrical coordinates 12 

0% = 27TOim P <Z 2 J[Ak] 2 k dk (3) 

^ 2e~ KZ °Zq sinh(fc d) . . 

k kKi ns cosh(k d ) + (kn v + 2 g Tj p/^)sinh(A: d) 

k v and Ki ns are the respective vacuum and insulator di¬ 
electric constants, while k is their average. And Qtf — 
2nq 2 j k Spuddle (E) is the Thomas-Fermi screening wave- 
vector which depends on the average density of states 
(Eq. 2). Ak is the potential solved from Poisson’s equa¬ 
tion which accounts for the distance of the impurities 
(z 0 ) inside the oxide, thickness of the oxide (d) and the 
screening length (1 Solved self-consistently (be¬ 

tween cfe and Fiddle) we determine the variance of the 
normal distribution of potentials (Fig. 2d). Over the 
dirty range, we can simplify it with a fitted equation 

&E ~ "b C (5) 


where C = 0.027eV 2 . This equation closely approximates 
the self-consistent calculation at the dirty limit. The vari¬ 
ation of cr m i n in presence of charged impurities allows us 
to quantify the competition between increasing modes 
and increased scattering. 

Analytical Model for cr m i n . The Landauer con¬ 
ductivity intuitively frames conduction as propor¬ 
tional to the transmission probability of electrons, T n , 
summed over all propagating and evanescent modes, 
where n is the mode index. cr min = GL/W = 


4 q 2 /h T n L/W. The general form for T n , deriv- 

n=0 

able by matching the pseudospinor wavefunctions across 

an n-p-n or p-n-p junction with barrier height Uo 
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yj(U o/Tive) 2 — q 2 and q n = nir/W is the transverse 
wave-vector in the channel that we sum over to get the 
total transmission. When k n is real then the trans¬ 
verse modes are propagating, while when k n is imaginary 
they become evanescent. Imaginary k n changes all the 
trigonometric functions to hyperbolic functions giving us 
an evanescent transmission T e = 1/ cosh 2 q n L when U Q is 
zero. An integral over a continuum of such cosh contribu¬ 
tions gives an overall factor of W/ttL which leads to the 
ballistic conductivity quantization (erg). For propagat¬ 
ing modes, the transmission probability T p picks up an 
additional scattering coefficient term from a series sum 
over the multiple scattering history, A/(A + L), where A is 
the electron mean free path in the presence of embedded 
impurities. The mean free path is vft sc where the mo¬ 
mentum scattering time r sc is determined from Fermi’s 
Golden Rule below. Combining all the elements in Eq. 
1, we arrive at the fan diagram in Fig. 1. 

Impurity scattering occurs through a 2D screened 
Coulomb energy, given at long wavelength by the 
Thomas Fermi equation, Vc(r) = Q 2 l(M^e^r)e ~ Kr . 
Using the pseudospin eigenstates, \Eqj(r) = 


l/y/2S (l e i9i ^y r e lki r normalized over area 
S', we calculate the scattering matrix element 
Vif = /d 2 r4>j(r)Uc(r)\Eq(r). In terms of scatter¬ 
ing wavevector and angle Ak = kf — ki, A 0 = Of — 0 $, 
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We can change to energy variables for elastic scatter¬ 
ing using |kf| = |ki| = E / (hv f) , (Ak) 2 = |kf — ki| 2 = 
k 2 + k 2 - 2kfki cos AO = 2E 2 (1 - cos AO)/(h 2 v 2 F ). For 
an impurity density ni mp and cross sectional area S (i.e., 
number of impurities nines'), Fermi’s Golden rule now 
gives us, h/r sc = |U/| di^E -£/&)( 1 cos Oj^)ni m pS. 

Converting sum into integral using the density of states 
(Eq. 2), and using the calculated expression for |U/| 2 
simplified for low energies, we get 
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The cosine integral followed by the delta function energy 
integral gives us 
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with Dp U ddie defined in Eq. 2. For E <C Acpft, 
the term in square brackets expands to 2 E 4 /h 2 v F K 2 + 
0(E 6 /h 4 v F K 4 ). We then get 


^ ^ Q ^imp-Spuddle /q\ 

t sc ~ 16 e ^ 2 

with k = ^ 2 D pudd i e /e 0 , giving us Ti/r sc = 
(^imp/lb)D pudd i e . Using the Einstein relation (diffusion 
coefficient V = v F r sc / 2), we get 
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At high impurity density, Dp uddle ~ 8a%/w 3 h 4 Vp (Eq. 
2). Using the approximate relation from Eq. 5 matching 
the self-consistent calculation fairly well in the dirty limit 
(Fig. 2), we get 
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Numerical model for cr m i n : We now show NEGF based 
numerical simulation results to calculate cr m i n in presence 
of charged impurities. We implement a discretized k.p 
Hamiltonian (H) to expedite computation. We use a 
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sequence of Gaussian potential profiles for the impurity 
scattering centers, 


^imp 

U(r) = ^2lJ n ex.p(-\r-r n \ 2 /2C, 2 ) (12) 

n= 1 


specifying the strength of the impurity potential at 
atomic site r, with r n being the positions of the impurity 
atoms and £ the screening length (~3 nm). The am¬ 
plitudes U n are random numbers following a Gaussian 
distribution with a standard deviation of lOOmeV 8 . This 
standard deviation is to be differentiated from the stan¬ 
dard deviation in the density of states description (Eq. 
2), which is a lumped description for the entire sheet in¬ 
stead of individual impurities. The Gaussian profile (Eq. 
12) is used to prevent the potential from going to in¬ 
finity at the scattering centers (Thomas-Fermi) and such 
approach is widely employed in the literature 8,15 18 . With 
U added to H , we calculate a m i n as a function ni mp (Fig. 
lb) by calculating average conductance over ^800 ran¬ 
dom impurity configurations. In the ballistic limit, cr m i n 
varies linearly with L/W, but as the sample gets dirt¬ 
ier, the cr m i n becomes less dependent on L/W. At high 
impurity limit, cr m i n becomes weekly dependent on ni mp 
and saturates around 4 q 2 /h. In most experiments, the 
device length L is larger than width W and therfore see 
a decreasing trend for cr m i n vs. ni mp such as Ref. 5 . The 
evolution of cr m i n from 4g 2 /nh to ~ 4 q 2 /h and therefore 
the missing i r can only be seen for devices with W » L. 
The differences between the numerical and the analyti¬ 
cal approaches most likely originate from the lack of ad¬ 
equate samples. 

Absence of Fabry-Perot as a signature of Klein tun¬ 
neling: Due to non-uniform doping along the metal- 

graphene-metal captured in our model by the differential 
dopings, a Fabry-Perot cavity is formed. Such a cav¬ 
ity leads to quantum interference oscillations and con¬ 
ductance asymmetry (n-n-n vs. n-p-n doping), seen in 
Fig. 3 in the ballistic limit. Such oscillations have been 
seen experimentally at low temperature in 2DEGs 19 , but 
are conspicuously missing for the lowest mode in single 
layer graphene (SLG), as seen in Fig. 3 left column. In 
contrast, the higher modes in the same column show os¬ 
cillations, as do all the modes for bilayer graphene (BLG) 
seen in the right column. The lowest mode in single 
layer graphene has forward and reverse propagating E — k 
bands with opposite pseudospin indices (bonding vs anti¬ 
bonding combinations of dimer p z orbitals) that disallow 
any reflection at heterojunctions. The resulting Klein 
tunneling 20 makes the heteroj unctions completely trans¬ 
parent to the lowest propagating modes and eliminates 
any Fabry-Perot oscillations. The parabolic lowest bands 
of BLG have twice the winding number around the Fermi 
circle (angle 2 6ij in the pseudospin eigenstate 4Lj) and 
thus a common pseudospin index, leading to finite reflec¬ 
tion and Fabry-Perot oscillations. 

We thus expect distinct behaviors of cr m i n vs. L/W 


in single layer and bilayer graphene. For large L/W , 



FIG. 3. (Color online) NEGF calculation of total conduc¬ 
tance G of single layer graphene and bilayer graphene reveals 
the nature of Fabry-Perot oscillation for the lowest mode. a,d 
shows linear and parabolic E — K in single layer and bilayer 
graphene. The lowest mode in single layer does not show 
any oscillation (b) but the bilayer does (e). The variation of 
minimum conductance and conductivity for single layer shows 
saturating G m in at 2 q 2 /h (c), while for bilayer graphene the 
minimum conductance never saturates and produces oscilla¬ 
tion m both Gmin and Cmin (f). 


G m in for SLG approaches 2 q 2 /h eliminating all tunnel¬ 
ing modes from source to drain and cr m i n = GL/W in¬ 
creases linearly (Fig. 3(a-c), already demonstrated in 
experiment 3 . For BLG (Fig. 3d-f), the conductance os¬ 
cillation for the lowest mode is manifested in the length 
dependence as well, leading to an oscillation in both G m [ n 
and cr m i n . For small L/W , the cr m i n saturates to 4 q 2 /(nil) 
and 2 q 2 /h for SLG and BLG respectively 13,21 . Such non¬ 
trivial transport behavior near the Dirac point is a mea¬ 
surable signature of Klein tunnel and reflection. 

Conclusion: The composite phase plot of graphene’s 
minimum conductivity is presented within a unified 
Landauer-Fermi’s golden rule and NEGF transport 
model. We show a general convergence of cr m i n vs. impu¬ 
rity concentration along with a quasi-saturation at high 
impurity concentration to ~ 4 q 2 /h irrespective of device 
dimensions. For high aspect ratios the increase in den¬ 
sity of states due to charged impurities results in a log¬ 
arithmically increasing cr m ^ n from the ballistic limit. On 
the other hand for low aspect ratios, the scattering due 
to charged impurities dominates and results in a power 
law decrease in the a m i n . For clean samples with con¬ 
ductance quantization, gating the sample into its low¬ 
est mode reveals a striking absence of low-temperature 
Fabry-Perot oscillations at low temperatures for SLG but 
not BLG, providing a signature of Klein tunneling. 
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